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1.  INTRODUCTION 


The  ability  to  provide  good  estimates  of  statistical  functions  of  ocean  data  is  an 
important  area  of  research  within  the  Naval  community.  A new  approach  for  estimating 
one  of  these  statistical  functions,  the  spectral  density  function,  is  the  focus  of  work  report- 
ed here. 

In  this  report  we  compare  two  discrete  time  spectral  estimation  techniques  for  esti- 
mating the  spectral  density  function  of  a continuous  time  process.  This  comparison  consists 
of  theoretical  results  which  give  the  asymptotic  performance  of  the  two  spectral  estimation 
schemes  along  with  simulation  results  for  finite  sample  size  estimates.  The  two  spectral 
estimation  schemes  considered  are  periodic  spectral  estimation  and  Poisson  spectral  estima- 
tion. In  the  first  case,  samples  of  the  continuous  time  process  are  taken  at  discrete  instants 
of  time  which  are  equally  spaced.  The  Poisson  spectral  estimate  is  based  upon  samples  from 
the  process  which  are  taken  at  points  of  time  determined  by  a stationary  Poisson  point 
process. 

Theoretical  asymptotic  results  for  the  bias  and  variance  of  the  Poisson  spectral  esti- 
mate have  recently  been  developed  [1, 2|  and  are  very  similar  to  the  asymptotic  expressions 
for  the  bias  and  variance  of  the  periodic  spectral  estimate.  Since  only  the  asymptotic  results 
are  available  for  the  Poisson  spectral  estimate,  the  performance  of  the  estimator  for  a finite 
number  of  samples  is  not  known.  Hence,  an  evaluation  of  the  performance  of  the  Poisson 
spectral  estimator  fora  reasonable  nu/nberof  samples  is  an  important  first  step  in  the  appli- 
cation of  the  theoretical  results. 
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n.  THEORY 


In  this  section  we  present  the  asymptotic  behavior  of  the  bias  and  variance  of  the 
periodic  and  Poisson  spectral  estimation  techniques.  The  periodic  spectral  estimator  is  a 
consistent  estimator  of  the  aliased  version  of  the  true  spectral  density.  The  Poisson  spectral 
estimator  is,  however,  a consistent  estimator  of  the  true  spectral  density  (3,4).  Similarities  be- 
tween the  asymptotic  performance  of  the  two  estimators  are  noted.  Derivations  of  these 
results  are  presented  elsewhere  [ 1 ] . 

Throughout  this  report  we  let  X = {X(t),  - '»<t<oo}bea  continuous  parameter, 
zero  mean,  second  order,  stationary  random  process  with  correlation  function  C(t)  eLj  and 
spectral  density  function  0(X)  given  by 

OO 

0(X)  = j C(t)e-'^‘dt/2jr  . 

— OO 

A spectral  window,  used  to  smooth  the  data  used  in  the  spectral  estimate,  is  generated  as 
follows:  Let  H(X)  be  a real,  even  function  which  satisfies 


/ H(X) 


dX  = 1 


/ IH(X)| 


dX  <<» 


The  spectral  window,  Wj^(X),  is  defined  by 
Wn(X)  = H (M^X)  , 
where  > 0 is  such  that 

Mn  ^ °°  and  Mj^/N  ^ 0,  as  N ^ . 

The  corresponding  covariance  averaging  kernel,  W]y|(t).  is  given  by 
W]^J(t)  = hft/Mjsi) 


where  li(t)  is  the  Fourier  transform  of  H(X), 


h(t)  = J H(X)e' 
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A.  Periodic  Spectral  Estimation 


Let  the  process  X be  defined  as  X = | X(t,^)|'”  where  the  sample  instants  |t,^| 
are  given  by 

tp  = n/p  ; n = 0,  ±1 . ±2, . . . 

and  p > 0 is  the  sampling  rate.  It  then  follows  14)  that  X is  a wide-sense  stationary  process 
with  zero  mean,  covariance  sequence 

c„  = C (n/p)  , n = 0.  ±1. . . . 


and  spectral  density  function 


s(X)  “ ^ ^ 0(pX  + 27rkp)  . 


Note  that  s(X)  has  period  27r  and  is  integrable  over  [-tt.tt]  with  Fourier  coefficients 


e“’'^s(X)dX,n  = 0.  ±1 


Now  the  spectral  density,  s(X),  of  X can  be  estimated  by  the  usual  technique  of 
smoothing  the  periodogram  of  the  discrete-parameter  stochastic  process  X.  The  estimate  of 
s(X)  based  upon  the  samples  ]X(n/p)||^_j  is 


sn(^)  = j Wj,j(X-p)l[^(p)dp 


where  Wj^(X)  is  the  spectral  window  (11-1 ) and  Ij^j(X)  is  the  periodogram 


^ X(k/p)e-ikX 


We  note  that  aliasing  occurs  in  (11-4)  unless  (/>(X)  is  zero  for  X e |-jrp,7rpl , i.e..  (piX) 
cannot  be  uniquely  recovered  from  s(X)  unless  (piX)  = 0 outside  of  [-7rp,7rp] . For  the  spe- 
cial case  that  0(X)  is  zero  outside  of  |-7rp,7rpl . then  0(’)  can  be  recovered  from  (11-4)  and  is 
given  by 


0(X)  =-^s(X/p)  . 
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In  the  case  in  which  the  process  X is  bandlimited,  say  from  -W/2  to  W/2,  then  the  sampling 
rate  p can  be  chosen  large  enough  so  that  0(*)  is  zero  outside  of  l-7rp,  ptt]  , i.e.  up  > W/2  or 


p > W/2jr 


The  quantity  W/27r  (cycles  per  second)  is  called  the  Nyquist  rate  of  sampling  and  is  the 
minimum  rate  at  which  one  can  sample  and  still  recover  the  spectral  density  of  the  original 
process  X.  If  X is  not  bandlimited  or  the  sampling  occurs  at  a rate  below  the  Nyquist  rate, 
then  aliasing  occurs  and  (/>(X)  cannot  be  recovered  from  s(X). 

The  form  of  (11-7)  motivates  the  estimate  for  the  case  when  aliasing  occurs.  Specifi- 
cally, define 


l//(X)  = ^s(X/p)  = ^ (^(X  + 27rkp)  . (11-8) 

k=-oo 

We  note  that  !//(•)  is  periodic  with  period  2vp  and  integrable  over  [-7rp,7rp] . The 

N 

estimate  of  i/^(X)  based  on  the  observations  {X(n/p)}^_j  is  given  by 

=7"sn(X/p)  (11-9) 

and  by  (II-5) 

OO 

/ Wn  (^-1^)  lN(M)dp  . (11-10) 


Note  that  since  i//fyj(X),  the  estimate  for  vKX),  is  periodic  with  period  27rp  then  at  best  we  can 
recover  i//{X)  from  the  sampled  data.  However,  as  the  sampling  rate  p tends  to  infinity 
i//(X)  tends  to  0(X). 

The  form  of  the  estimate  (11-10)  for  vi'isj(X)  is  convenient  for  obtaining  asymptotic- 
properties  of  the  estimator,  however,  for  practical  implementation,  the  following  form  is 
used 


/ 

1 


1 

27r  Np 


I 


j=l 


N-1 


2 

m=l 


N-m 


(11-11) 

where  h is  defined  by  (11-3).  The  above  form  of  \l/'^  is  obtained  by  substituting  (11-6)  into 
(II-IO)  and  making  use  of  (Il-l ) and  (11-2). 


-I 

1 


'i 

i 


We  now  look  at  the  properties  of  the  estimate  The  effects  of  aliasing  will  be 

emphasized. 

Theorem  1 . Let  the  covariance  function,  C(t),  of  the  process  X satisfy 

(i)  tC(t)eL, 

OO 

(ii)  ^ n |C(n/p)|  <«>  . 
n=l 

Then  the  asymptotic  mean  of  the  estimate  is  given  by 

OO 

E (•^)]  (//(Pldp  + O (-jjj-) 

~oo 

as  N OO  and  the  0(1  /N)  term  is  uniform  in  X. 

The  proof  of  Theorem  1 is  given  in  ( I J . By  letting  N oo  we  see  that 

lim  E r^'isjfX)]  = i//(X) 

N ->  OO  L J 

and  the  estimator  i/'{yi(X)  is  an  asymptotically  unbiased  estimate  of  ipCX).  However,  recall 
that  we  are  interested  in  estimating  0(X),  Ihe  spectral  density  of  X.  and  not  \p{\).  the  spec- 
tral density  of  X.  From  (II-7)  and  (11-9),  the  best  that  we  can  hope  for  is  an  estimate  of 
^(X),  which  contains  an  aliasing  error,  e(X).  defined  by 

€(X)  = i//(X)-0(X) 

OO 

= ^ 0(X  + Ztt  kp)  . 
k=-oo 
kstO 

Note  that  the  aliasing  error  is  independent  of  N,  the  total  number  of  samples  used  to  form 
the  estimate,  and  therefore  cannot  be  reduced  by  increasing  the  sample  size. 

The  following  theorem  gives  the  covariance  of  the  estimate  i^fyj(X)  for  Gaussian 
processes. 
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Theorem  2.  Let  X be  a Gaussian  process  with  covariance  function  C(t)  satisfying 


(i)  tC(t)eL, 

OO 

(ii)  ^ n|C(n/p)|<< 


n=l 

Let  the  spectral  window  be  given  by  (H-l)  and  suppose  H(X)  is  bounded  and  differentiable 
with  bounded  derivative.  Then  asymptotically  (as 


cov 


M 

+ j H^fiOdu 


'M 


N. 


M 


N 


where  ^ is  the  27rp-periodic  extension  of  the  Kronecker  delta  function 


lo, 


1 , if  X-/i  (mod  2vp) 
otherwise  , 


and  the  O(')  term  is  uniform  in  X and  p.  Moreover, 

M 


Var 


2ir-~  [(t>(K)  + e(X)]^  j H-(u)du  , X ¥=  0 (mod  Trp) 


Mn  ~ r ^ 

47r^^  [(^»(X)  + e(X)]  ^ J H-(u)du  , X = 0 (mod  7rp) 


The  proof  of  this  theorem  follows  the  proof  in  [5,  Theorem  5.6.2] . 

The  results  of  Theorem  2 are  valid  for  non-Gaussian  processes  as  well,  under  addi- 
tional conditions  on  the  fourth-order  cum  ulant  function  of  the  process  X [5,  6J.  It  can 

/V 

now  be  seen  by  Theorems  1 and  2 that  for  Gaussian  processes  the  estimate  i//j>^j(X)  con- 
verges in  the  mean-square  sense  to  i/^fX)  as  N^oo.  However,  it  does  not  converge  to  0(X) 
due  to  the  aliasing  error  e(X). 

B.  Poisson  Sampling 

The  motivation  for  Poisson  sampling  is  discussed  in  [1,2].  Among  the  sampling 

N 

.schemes  which  yield  consistent  estimates  of  0(X)  from  the  observation  {2((t[^))|^_l . Ilic 
Poisson  sampling  scheme  is  the  simplest  and  the  most  fundamental  1 1 | . Here  we  only  intrcv 
dQce  the  estimate  and  state  its  bias  and  covariance  for  the  sake  of  comparison  with  the 

/N 

periodic  estimate  >//[«^(X). 
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Assume  that  tlie  process  X is  sampled  at  instants  {tj^ } given  by 


1 

i. 

i 

t 

[ 

f- 


I 


{ 


tn  - 0 


tn  = ‘n-1  +“n  ■ »= 


(11-12) 


where  the  {Oj^jare  independent  identically  distributed  positive  random  variables  with  a 
common  exponential  distribution  F(x)  = 1 It  is  assumed  that  the  sampling  instants 

OO 

independent  of  the  process  X.  Note  that  /3  is  the  average  sampling  rate  and  as 
shown  in  [2],  almost  every  realization  of  the  point  process  {tp}~_Q  has  a “density”  of 
points  /3. 

N 

Given  the  observation  {X(t|^)}j^_j,  we  consider  an  estimate  (j>^{X)  ot  the  spectral 
density  0(X)  as  follows  [ 1 1 : 

N-1  N-n  j 

y X(tk+„)X(tk)h/-^^^'\cosX(t,,  + „-t^)  (11-13) 

n^l  k=l  \ / 

where  h(t)  and  are  as  in  (Il-l ) and  (11-3).  It  can  be  seen  that  the  estimate  has  an  almost 
identical  structure  to  that  of  the  periodic  estimate  (11-1 1 ).  Define  the  bias  of  the  estimate 
(11-13)  by 

b = E pN^^)]  • 

We  then  have  the  following  results  1 1 1 stated  here  without  proof. 

Theorem  1':  Let  the  covariance  function  of  the  process  X satisfy  tC(t)  e L| . Then, 
asymptotically  as  N^o®, 

OO 

E ~ / Wf^(X-/a)  0(ju)d/i  + 0 

— OO 

where  the  0(  1/N)  term  is  uniform  in  X. 


Theorem  2':  Let  X be  Gaussian  with  covariance  function  C{t)  satisfying  |C(t)|  < 
b(t)  for  all  t,  where  b(t)  denotes  any  continuous  even  nonnegative  function  on  the  real  line 
which  is  nonincreasing  over  [0,°°)  such  that  tb(t)e  L| . Let  the  spectral  window  be  given  by 
(Il-l ) and  suppose  H(X)  is  bounded  and  differentiable  with  bounded  derivative.  Then 
asymptotically  as  N-»-oo 


2 

cov  j^0j^(X),  0j^(Ai)]  = 27r/3(6j^  ^ + 6;^  _^)  j^0(\)  + ^^j  j H2(u)du 

-iOO 

/“n  Mn\ 


where  ^ is  the  Kronecker  delta  function 
1 , X=/i 

®X,u  “ 

0,X^p  , 


and  the  0(‘)  term  is  uniform  in  X and  /i.  Moreover 


M 2 ^ 

Var  p^fX)]  [*(X)  + ^]  (i+So,X>  / H2(u)du  . 

~oo 

The  results  of  Theorem  2'  are  valid  for  non-Gaiissian  processes  as  well,  under  addi- 
tional conditions  on  the  fourth-order  cumulant  function  of  the  process  X 1 1 1 . 

C.  Comparison 

In  order  to  compare  the  asymptotic  performance  of  the  periodic-estimate  (11-1 1 ) 
and  the  Poisson-estimate  (11-13),  we  assume  the  two  estimates  to  have  the  same  number  of 
samples  N,  identical  window  functions  H(X)  and  identical  sampling  rates  p=(3.  It  is  then 
seen  by  Theorems  1'  and  2'  that  the  Poisson-estimate  0[y[(X)  converges  in  the  mean-square 
sense,  as  N-«o,  to  the  spectral  density  0(X)  for  any  positive  value  of  the  average  sampling  rate 
0.  On  the  other  hand,  it  is  seen  by  Theorems  1 and  2 that  the  periodic  estimate  (11-1  1 ) 
does  not  converge  to  0(X),  as  unless  0(X)  = 0 for  |X|  > np.  Comparing  now  Theorems 
1 and  1 , we  find  that  for  the  periodic-estimate  0|q(X)  we  have 

OO  OO 

E [iN(X))  = / (^)J  0(M)d/i+  J |^-iwN(-^)]  e{p)dp 


where  the  second  integral  represents  the  increase  in  the  mean  of  the  estimate  due  to  the 
aliasing  error.  On  the  other  hand,  for  the  Poisson  estimate  0jq(X)  we  have 


111.  GENERATION  OF  DATA 


In  this  report  we  consider  three  spectral  density  functions  which  are  associated  with 
different  second  order  random  processes.  Specifically,  the  spectral  density  functions  consid- 
ered are 


«“  + x~ 


^>(X)  = 


+ 4a^X^ 


(Ill-la) 


(lll-lb) 


0(X)  = 


2 aX^/ir 


2 aX^/TT 


+4d^y 


(x2-x,2) 


2 . T..  2 

4a"X 


(III-lc) 


These  spectral  densities  have  been  chosen  because  they  possess  certain  characteristics  which 
are  found  in  data  obtained  from  the  real  world.  The  first  spectral  density  (Ill-la)  is  a first 
order  Gauss-Markov  process.  This  spectrum  is  fairly  fiat  and  could  represent  wide-hand  noise 
which  has  a constant  spectrum.  The  second  spectral  density  function  (lll-lb)  could  repre- 
sent a carrier  modulated  by  a narrowband  signal.  Finally,  the  last  spectral  density  function 
(III-lc)  represents  two  closely  spaced  carriers  modulated  by  narrowband  signals. 

The  remainder  of  this  section  describes  the  generation  of  data  which  represent 
sampled  values  of  the  associated  random  processes.  Tlie  aliased  version  of  each  of  the  three 
spectral  densities  is  derived  and  the  dynamic  representation  of  the  three  spectral  densities  is 
obtained.  The  dynamic  representation  provides  a recursive  method  of  generating  the  sample 
values.  Finally,  the  statistics  of  the  random  variables  which  are  needed  in  the  dynamic- 
representation  are  given. 

A.  First  Order  Gauss-Markov  Spectral  Density 


The  first  order  Gauss-Markov  spectral  density  function  (lll-la)  has  the  simplest 
structure  (Figure  1 ) of  the  three  spectral  densities  used  in  the  simulation.  It  represents  a 
nonbandlimited  low  pass  spectrum  and  provides  a test  of  the  bias  of  the  spectral  estimators 
(11-10)  and  (II-l  2).  The  correlation  function  R(t),  is  given  by 

R(t)  = . 

Since  the  spectral  density  function  is  not  bandlimited  and  the  mean-square  bandwidth  also 
does  not  exist,  we  will  refer  to  the  half  power  point  (or  3 dB  down  from  the  maximum)  as 


where  f|>(t-T)  is  given  by  (A-5). 

Therefore,  a recursive  representation  is  obtained  for  the  process  by  replacing  t with  t|^^.j  and  ^ 

T by  t|^  in  equation  ( III-2).  This  yields  | 


1 


/ 

i 


with 


(v)  di^ 


In  order  to  simulate  the  X process  thef  t|^ } must  be  known  and  the  ej^  must  be  ob- 
tained. Let  XCtg)  = X(0)  = 0.  Tlie  ej^  are  zero  mean  Gaussian  random  variables  since  they 
are  obtained  by  a linear  operation  (integration)  on  Gaussian  random  variables  u(t).  Now 
only  the  variance  of  e|^  has  to  be  obtained  in  order  to  completely  specify  the  behavior  of  the 
X process  at  the  instants  {tj^ }. 

The  variance  of  is  computed  as  follows 


*k+l  /k+1  _ 


= I f 


k+1  ^^k+1 


= / / 
*k 


e ‘ e ‘ u(t/)Li(r})di'dTi 


-Qfltk+i-J')  -“(tk  + i-i?) ^ 
e ' e ‘ 8(i^-rj)di’dv 


.‘k+1 


-/  '■ 


^-:a(tk+|-i') 


di' 


The  variance  of  = G C|^  is 


Etj^  = G-  • Ee“  = R(0) 


2ar(tk+i  - 


‘k^j 


(I1I-4) 


Also,  it  can  be  shown  that  the  cj^  are  independent  random  variables.  This  fol- 
lows from  the  definition  of  e|^  in  terms  of  the  white  Gaussian  noise  which  makes  the  C|^ 
uncorrelated,  and  since  the  e|^  are  Gaussian  it  follows  that  they  are,  in  fact,  independent. 

Poisson  Sampling  Model 

The  following  model  de.scribes  the  method  used  to  generate  the  sampling  points  and 
the  values  of  the  process  at  the  sampling  points.  Let  the  sampling  instants  {tj, } be  generat- 
ed by  a Poisson  point  process  (11-12). 
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Then  the  dynamic  representation  is 


X(0)  = 0 

X(tk+i)  = e ‘‘^'xUki  + Tik 

where  the  t^  are  given  by  (II-I2)  and  “ sequence  of  independent  zero  mean  Gaussian 
random  variables  with  variance  (equation  (III-4)) 

. r|.e--<'k+|-‘k>] 

% L J 

Note  that  the  variance  will  change  for  each  k. 

Periodic  Sampling  Model 

In  the  model  for  the  periodic  sampling  case  the  {t^ } are  determined  by 

tk  = kh.  k = 0,  1.  2 (I1I.5) 

where  h is  tlie  time  between  samples.  The  sampling  rate,  p.  is  given  by 
P = 1/h  . 

Let  Xk  = X(kh).  Then  the  dynamic  representation  simplifies  to 
Xo  = 0 

Xk+,  = e-«''X^  + 0 


k^'^k 


where  the  {^k^  independent  identically  distributed  zero  mean  Gaussian  random  varia- 
bles with  variance 

aj  = [l-e-2/P|  . 

Note  that  the  variance  of  ^k  does  not  depend  upon  k. 

B.  The  Narrowband  Spectral  Density 

The  spectral  density  used  to  represent  the  process  X for  the  narrowband  case  is  given 
by 


0(M  =~ 


+ 4a^X 


2v2 


(111-6) 
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where  2 a/rr  is  a normalization  constant  such  that  the  correlation  I'unction  evaluated  at  0 is 
equal  to  1 . i.e.,  R(0)  = 1 . A representative  form  of  the  spectral  density  given  by  equation 
(111-6)  is  shown  in  Figure  2.  The  correlation  function  is  given  by 


= f 0(X)e‘‘^dX 

= e~^*  [cos  Wj,t sin  w^jt).  t > 0 
\ '^o  ' 


(111-7) 


where 


2 ,2  -) 
w = X - a~ 
o o 


The  maximum  of  the  spectral  density  occurs  at  X = X^  and  is 
= l/2jra  . 

The  half-power  point,  X|,,  is  obtained  by  solving 
0(X,,)  = 0(Xo)/2 
for  X|,  which  yields 


(111-8) 


(111-9) 


By  adjusting  Xq  and  X|,,  the  “bandwidth”  of  the  narrowband  stochastic  signal  can  be  made 
to  provide  an  appropriately  narrow  peak  to  the  spectral  density. 

The  aliased  version  s(X)  of  (111-6)  is  derived  as  follows.  Let  r(n)  be  the  correlation 
function  of  the  discrete  random  process  X.  Then  r(n)  = R(n/p)  where  R(')  is  the  correla- 
tion function  of  the  continuous  time  random  process.  Therefore.  s(X)  is  given  by  the  Poisson 
summation  formula  as 


s(X)  = ^ R(n/p)e  • 


(111-10) 


Since  R(‘)  is  an  even  function,  s(*)  may  be  rewritten  as 


s(X)  = ~ R(0) 


+ ^ R(n/p)e~‘”^+  ^ R(n/p)e‘”^ 


= R(0)  + 2Re  ^ R(n/p)e"‘"^ 


From  (111-7)  the  correlation  function  can  be  written  as 

, Wo+ia  iT(w,,+ia)  _ w^-ia  irf-w^^+ia) 
-w„  2w^ 


2 £ ^inl(w„+l.,/.-XI 

11=1  ° 11=1 

'^o“'“  V in[(-w„-Ha)/p-Xl 

Z 


2wo 

n=l 

Wo+ia 

e(iWo-a 

2wo 

g-a/p 

'’w 

“”o 

( Wo+ia) 

i/p-iX  WQ-ia  (iwQ-a)/p-iX 
- ^ ^ ^ 

i)/p-iX  2w_,  , (-iw,,-a)/p-i 

o 1 _e  ^ 

cos  B~-t-i  sin  B~-e~‘‘^^ 

1 +e~~''^lP  -2e~'^^P  cos  B“ 


. cos  B'*^-i  sin 

<'^«'‘"^+e-2a/P-2e-"/^cos  B+ 


where 


B"^  = w^/p  + X and  B"  = w^/p  - X 
Therefore. 


Re  ^ R(n/p)e~'”^ 

11=1 

cos  B“-e~“^^  sin  B~  cos  B^-e“^^^  - — sin  B"*^ 

= j.-a/p  2 + 2 

1 +c~-“^^-2e~^/^  cos  B"  1 -He~2a/P_2c"^/^  cos  B^ 


Hence,  the  aliased  spectral  density  is 
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‘ Define  the  random  variables  and  i7|^  as 

M ‘k+l 

[ ^ j 011  (tk+i-T)u(T)dT  , k = 0,1,2,... 

r 

! ‘k+l 

[ T?jj  = y 021  (‘k+l“^)  ^ ^ 0,1,2 

^k 

and  assume  the  initial  values  of  X and  X2  are  zero,  i.e., 

X(0)  = X2(0)  = 0 . 

[ 

i Since  u(t)  is  a zero  mean  white  Gaussian  noise  the  random  variables  and  are  also  Cans-  \ 

* 00  00  *• 

r sian  random  variables  with  0 means.  The  sequences  {di, } and{Tji,}  are  each  a se- 

j k=l  k=l 

quence  of  independent  random  variables  since  each  (and  tj|^)  is  generated  by  white  Gaus-  ] 

I sian  noise  over  distinct  time  intervals.  However,  is  correlated  with  since  the  time  I 

^ interval  of  noise  which  generates  each  random  variable  is  the  same.  Note  that  dj,  and  for  j 

t k=^n  are  independent  by  the  same  reasoning  as  above.  ■ 

! In  order  to  completely  describe  d^.  and  1?]^  an  expression  for  their  variances  must  be  i 

computed.  The  spectral  density  for  the  white  Gaussian  noise.  u(t),  is  constant  and  equal  to  j 

1 

>'  I 

I 0u(X)  = 2a/7r  j 

I 

I and  the  correlation  function  of  u(t)  is  given  by  j 

I ^u^^^  ~ ‘ 2a6(t)/7r  = 4a6(i/  . j 

The  variance  of  d|^  is 


Similarly,  the  variance  of  is 


2 ^2  ^2  -2a Ac  / 

a„_  = X - X _ e ^ h 


o o 


^awQ  sin  2 A|^  -a“  cos  2 Aj^  + X^  )/'''q  (III-20) 


The  cross  correlation  of  6^  and  is  given  by 


ri^  = E i7|^  = 4a  j 0j  j (t)  (^21  ‘It 


9 -2aAi,  2 

= aX“e  (cos  2 Wq  A|^-l)/w^  , 

and  the  correlation  coefficient  is 


(III-21) 


We  now  have  the  dynamic  representation  for  the  narrowband  process  X(t),  with 
initial  conditions  t^  = 0,  X(0)  = X2(0)  = 0 given  recursively  by 


X( t|^^ I ) - 0 1 1 ( A|^)  X ( t j^)  + 0 j T ( Aj^)  X')(t]^)  + d|^ 
X-i(t|^^  (A|^)X(t|^)  + ^'i')  ( A|^)  1 


(III-22ai 


(lll-22b) 


where  A|^  = t|^+j  - t^.,  the  are  given  by  equations  (III-17a)  --  (III-17d).  and  the  and  tJ|^ 

■*>  1 

are  correlated  zero  mean  Gaussian  random  variables  with  variances  Oa  - and  a_  “ given  by 

''k 

(111-19)  and  (III-20)  with  correlation  coefficient  p|^  (111-21). 

The  generation  of  the  narrowband  data  is  begun  by  transforming  two  sequences  of  in- 
dependent zero  mean  Gaussian  random  variables  into  two  sequences  of  correlated  random 
variables  with  a specified  cross  correlation.  This  transformation  will  provide  the  proper  meth- 
od of  generating  sequences  of  the  random  variables  tj^.  and  0|^.  Given  two  sequences  of  zero 
mean,  unit  variance  Gaussian  random  variables  {X|^}and  {Y|^}  we  will  obtain  the  transforma- 
tion which  will  yield  two  sequences  of  Gaussian  random  variables  } and  with  vari- 
ances ’ '■espectively,  and  correlation  coefficients  {pp;}. 


For  simplicity  of  analysis  the  time  subscript,  k,  will  be  omitted.  Let 


W "cos  a - sin  a"! 

.Zj  sin  a cos  a J I t?  1 


(111-23) 


where  0,  j?  are  the  desired  Gaussian  random  variables.  Note  that 


0~  'cos  a sin  a 1 'W' 

_-sin  a cosaj  Z_ 


(111-24) 


We  will  find  the  appropriate  a = a(ai , (77,p)such  that  W and  Z are  independent.  Also,  we 
2 2 

find  the  variances,  Oy^  and  of  W and  Z. 


(i)  Obtain  from  a zero  mean,  unit  variance  Gaussian  random  number  generator  two 

sequences  of  numbers  and  Also,  obtain  a realization  of  the  sampling  | 

points  {t|^}^_|.  For  the  data  used  in  the  periodic  spectral  estimate  the  set  {t|^^}  is  obtained  | 

by  (111-5)  and  for  the  Poisson  case  by  (11-12).  j 

(ii)  Compute  Var  . Var  , and  from  (111-19),  (111-20)  and  (111-21)  | 

respectively. 

(iii)  Compute  Var  I Wj^l , Var  (Z^| , and  from  (111-27),  (111-26)  and  (111-25).  | 

(iv)  Set  W^=  fVarlW^I  i ' X,^  and  Z,^  = {Var  | Z,^! } ' ^ 

(v)  Compute  6/ and  j 

(vi)  F'inally,  compute  X-i(tj^^l ) from  (lll-22b)  and  then  the  observation  X(t|^+I ) by  j 

(lll-22a).  j 

' i 

i 


* 


22 


C.  TWO  CLOSELY  SPACED  NARROWBAND  SIGNALS 


% 


[ 


IF 


A spectral  density  resembling  two  closely  spaced  narrowband  stochastic  signals  is  the 
final  spectral  density  case  investigated  (see  Figure  3).  This  data  set  demonstrates  the  ability 
of  the  Poisson  spectral  estimation  algorithm  to  resolve  two  narrowband  signals  closely  spaced 
in  frequency.  The  spectral  density  is  of  the  form. 


<>(X)  - 


2 aX“/jr 


2 aX-/7r 


+ (x2-X|2)2  t4A2 


(in.28> 


where  Xq  and  X|  are  the  center  frequencies  of  each  narrowband  signal  and  a is  as  previously 
defined. 

In  this  case  we  assume  that  Xq  and  X|  are  two  zero  mean  independent  random 
processes  with  correlation  functions 


1 


I 

i 


(T)  = e“^^  (cos  W:  T - — sin  W;  T ) ,i  = 0,l 

Xi  ^ 1 Wj  1 / 

where 

2.2  -> 

Wj  = Xj  -a-  i = 0,1 

and  corresponding  spectral  densities 


OO 

<|>j(X)=  J"  e'^’’R^.  (r)  dr  = 


2 a X~ln 


(x2-x/)-.4a^X^ 


j = 0.1 


Define  the  process  X(t)  as 

X(t)  = Xo(t)  + X|(t). 


The  aliased  version  of  <t>(X)  is 


i//(X)  = ^0^^^  + '/'l(X) 
where  i//q(X)  and  i//](X)  are  as  in  (Ul-1 1)  . 


<//i(X) 


A r 1 W:  a O wn  \ ->  2W: 

l-q^+2  (q  -q)cos— (q'^+q)  sin  * cosi^+ 2 q"^  sin — i. 

^ _L  L P w- pJ  P Wj ^ 

27rp  W:  / 2w.  2X\ 

(l+q'^)  - dfq-’-q)  cos cosA+ 2q“  (cos  — i+ cos — I 


where  q = eT^/P,  i=0, 1 


All  other  parameters  are  as  previously  defined. 


The  Dynamic  Representation 

Since  the  spectral  density  ^(X)  is  of  the  form 


Po(iX) 

2 

+ 

P((iX) 

Oo(iX) 

Q](iX) 

then  the  dynamic  representation  of  the  process  X(t)  is  simply  the  addition  of  the  dynamic 
representations  for  Xgft)  and  X|(t).  These  dynamic  representations  are  derived  in  exactly 
^he  same  manner  as  for  the  narrowband  case  considered  in  the  last  section. 

The  process  X(t)  for  this  case  is  generated  from  (111-22)  . 

X(tk+i)  = Xo(tk+i)  + X|  (tk+i) 

where 


Xi  (tk+i)-  0i,ii  ^*k>'^<^i,12^^k^  ^i,2^^k^'^^i,k 

i=0,l 

^i,2  <‘k+l  > *^i.21  ^^k^  ^i  ^*k>  K22  ^^k^  ^i.2  ^‘k^  ’’i.k 
In  order  for  Xq  and  X|  to  be  independent 

E l®0,k  ^l.kl  ""  ^ I^O.k’^l.k*  = E ieo,k77i,kl  = E t’7o,k®l,k'  " ® 

must  be  satisfied  and  hence  four  sequences  of  random  numbers  must  be  generated  instead  of 
two  sequences  as  before.  The  matrices  jdt)  and  0]  jj(t)  for  this  case  are  given  by  (111-17) 
where  Wq  = Xq"^  - a*"  for  0q  jj(t)  and  w]"  = Xj-  - a“  for  (>  j jjd). 


The  expressions  for  6^  and  tjj  are  taken  from  (111-24)  and  are  of  the  form 
®i,k  = \k  “i,k  + Kk  sin  «i,k 


i=0,l 


^i,k  = -Ai,k  sin«i,k  + Bi,k‘-‘«sai,k 


where  {Aj  and  {Bj  are  four  sequences  of  independent,  zero  mean,  unit  variance 
Gaussian  random  numbers  and  Oj  k ‘s  defined  by  (111-25). 


“i,k  = ‘an 


[2pia0.a^./ 


The  values  of  pj,  oq.  and  o^.  are  obtained  from  (111-19)  through  (111-21 ).  Specifically, 


2 = , + -2a  ^k  / „ 


Oa  = 1 + e 


1 'y  \ 1 

a Wj  sin  2 Wj  Aj.  + a“  cos  2 Wj  A|^  - Xj“  I /wj" 


^i 


-1  T -)  —2a  Ai,  / T 

“ “ X;  - Xj“  e I a Wj  sin  2 Wj  Aj.  - a“  cos  2 Wj  A^,  + X 


i")/'^i“ 


3j=  (aXj^e  “ 


-2a  A 


k 


(cos  2 Wj  A|^  - 1 )/ w 


i‘)^"®i,k%,k 


where 


Wj^  = Xj2  _ a2  j 
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IV.  SIMULATION  RESULTS 

The  simulation  was  performed  for  each  of  the  three  spectral  densities  using  the  uni- 
formly sampled  periodic  estimator  and  the  Poisson  estimator.  For  each  case  simulated  both 
estimators  were  provided  with  the  same  or  corresponding  parameter  values.  The  number  of 
time  samples  was  held  constant  and  equal  to  1000,  and  the  Parzen  window  h(t),  given  by 

l-6t2  + 6|t|^  It]  <7 

2(1 -Itl)'*  CV-l) 

0 KItl 

was  used  as  the  covariance  averaging  kernel  for  both  estimators.  The  Parzen  spectral  window 
was  chosen  because  it  reduces  the  variance  of  the  estimate  and  is  everywhere  positive  in  the 
frequency  domain.  This  spectral  window  eliminates  the  possibility  of  yielding  negative  esti- 
mates for  the  periodic  estimator  and  reduces  the  probability  of  negative  estimates  for  the 
Poisson  estimator.  Corresponding  sampling  rates  p for  the  periodic  case  and  (average)  i 

sampling  rates  for  the  Poisson  were  used  for  each  spectral  density,  specifically  p = |3  for  each  ] 

case.  Several  different  values  of  the  sampling  rates  were  considered  for  each  spectral  density 
in  order  to  determine  their  effect  on  the  estimates. 

The  value  of  the  truncation  point  M]q  of  the  estimators  was  chosen  for  each  case  as 
the  value  which  minimizes  the  mean  squared  error  over  the  frequency  range  investigated. 

The  procedure  for  choosing  the  was  to  compute  several  estimates  of  the  same  spectral 
density  while  varying  and  then  selecting  the  for  the  case  of  smallest  mean  squared 
error.  The  values  of  Mj^j  for  the  periodic  and  the  Poisson  estimators  do  not.  in  general, 
coincide.  This  is  because  the  periodic  estimator  is  actually  estimating  the  aliased  spectral 
density  rather  than  the  true  spectral  density  as  in  the  Pois.son  case. 

The  results  of  the  simulation  for  each  of  the  three  spectral  densities  are  presented  in 
the  following  sections  and  a summary  of  this  work  appears  in  [ 7 ] . 

A.  First  Order  Gauss-Markov  Case 

For  this  case  three  sampling  rates  were  used  for  both  the  Poisson  and  periodic  esti- 
mates. These  are  p = /3  = 2,  4 and  27r  radians/second.  For  the  Poisson  case  the  sampling  rates 
correspond  to  average  values.  These  rates  were  chosen  as  being  multiples  (2,  4,  2;r)  of  the 
half  power  point  of  the  spectral  density.  They  are  representative  of  sampling  at  approxi- 
mately the  4,  d and  1 3 dB  down  points  of  the  spectral  density.  The  number  of  data  points 
used  in  the  simulation  was  N = 1000.  Mj>j  was  selected  on  the  basis  of  smallest  mean  square 
error,  the  Parzen  spectral  window  (IV-2)  was  used  and  a=  1.0. 

Spectral  estimates  for  these  sampling  rates  are  shown  in  Figures  4,  5,  and  6. 

The  continuous  curve  appearing  on  the  graphs  represents  the  true  continuous  time  spectral 
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density  0(X).  The  spectral  estimate  ${\)  for  the  Poisson  case  and  ^(X)  for  the  periodic  is 
represented  by  a continuous  curve  with  square  symbols.  The  graphs  for  the  periodic  case  also 
contain  a dotted  curve  representing  the  aliased  spectral  density  i//(X). 

Figure  4 compares  the  two  sampling  techniques  for  p = j3  = 2 radians/second.  The 
periodic  graph  shows  a large  amount  of  aliasing  with  period  27rp  which  inhibits  for  any  X a 
reasonable  representation  of  the  true  spectral  density  <^(X).  On  the  other  hand  the  graph 
representing  the  Poisson  technique  does  not  indicate  any  aliasing  effects  and  for  all  X pre- 
sented on  the  graph  a fair  estimate  is  provided. 

By  doubling  the  sampling  rate  to  p = /3  = 4 it  is  seen.  Figure  5,  that  aliasing  in  the 
periodic  estimate  is  reduced.  This  is  because  as  p becomes  larger  0(X)  becomes  smaller  out- 
side the  interval  [-rrp,  rrp] . This  estimate  for  the  periodic  case  is  still  very  poor.  The  corre- 
sponding Poisson  spectral  estimate  is  very  good  with  only  slight  fluctuations  from  the  true 
spectral  density  </)(X).  The  reason  this  estimate  is  so  good,  especially  at  the  origin  (X  = 0).  is 
that  from  (11-14)  ^ = 4 rad/s  is  near  the  optimal  (rr  rad/s)  forX  on  or  near  the  origin. 

The  periodic  spectral  estimate  for  the  sampling  rate  p = |3  = 27r  rad/s,  Figure  6, 
begins  to  resemble  the  true  spectral  density  </>(X)  but  is  still  heavily  influenced  by  the  aliasing 
error.  The  Poisson  estimate  still  provides  a good  estimate,  especially  in  the  tail 
X>  2. 

For  curiosity’s  sake  the  frequency  axis  of  the  Poisson  graphs  was  extended  to  X = 21 
and  spectral  estimates  were  computed.  The  resulting  estimate,  which  is  not  shown  here, 
closely  followed  the  true  spectral  density  all  the  way  to  X = 21  with  no  sign  of  aliasing  or 
increasing  error. 

B.  Narrowband  Case 

Here  the  two  estimates  (II-l  I)  and  (11-13)  were  evaluated  for  the  sampling  rates 


p = 13=  lira  = 0.6  rad/s 
p = |3  = 2.2  rad/s 
p = d = 3.2  rad/s 


(IV-2a) 

(IV-2b) 

(IV-2c) 


The  sampling  rates  were  chosen  to  represent  a sampling  rate  well  below  Nyquist  (IV-2a),  one 
at  the  3 dB  point  on  the  spectral  density  (IV-2b)  and  finally  a sampling  rate  (IV-2c)  far  above 
the  center  frequency.  The  first  rate  (27ra  rad/s),  in  addition  to  being  below  Nyquist,  is  well 
below  the  peak  of  the  spectral  density.  From  (11-14),  this  is  the  optimal  rate  for  estimating 
the  peak  at  X = Xq. 

As  with  the  first  spectral  density  simulated,  N = 1 000.  was  selected  on  the  basis 
of  smallest  mean  square  error,  the  Parzen  spectral  window  (IV-2)  was  used,  and  the  parameters 
a and  Xq  were  chosen  to  be  2.1/22  and  1.0  respectively. 
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0.0  1.0  2.0  3.0  4.0  5.0  6.0 

FREQUENCY  <rad/s) 

Figure  4a.  Poisson  spectral  estimate  of  first  order  Gauss-Markov  process 
sampled  at  (3  = 2 rad/s  with  Mjsj  = 7. 


0.0  1.0  2.0  3.0  4.0  5.0  6.0 

FREQUENCY  (rad/s) 

Figure  4b.  Periodic  spectral  estimate  of  first  order  Gauss-Markov  process 
sampled  at  p = 2 rad/s  with  = .S. 
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FREQUENCY  lrad/s| 

Figure  6a.  Poissrui  spectral  estimate  of  first  order  Gauss-Markov  process 
sampled  at  d = 27r  rad/s  with  = 6. 
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Figure  6b.  Periodic  spectral  estimate  of  first  order  Gauss-Markov  process 
sampled  at  p = 2fr  rad/s  with  Mj^  =15. 
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The  results  for  the  simulation  of  this  spectral  density  are  presented  in  Figures  7 
’ through  9.  Figure  7 shows  the  results  of  the  estimates  for  p = 0 = 27ra  rad/s.  This  sam- 

pling rate  is  well  below  the  band  which  is  centered  at  = 1 . The  periodic  estimate  shows 
« no  resemblance  whatsoever  to  the  true  spectral  density  (^(X).  The  majority  of  the  error  is  due 

to  spectral  aliasing.  The  Poisson  estimate  for  this  case  provides  a slight  bias,  in  frequency,  but 
otherwise  a very  good  estimate  of  the  spectral  density.  No  aliasing  is  present  in  this  graph  for 
any  X.  Figure  8 shows  the  spectral  estimates  for  the  second  sample  rate  p = /J  = 2.2  rad/s. 

The  Poisson  estimate  in  this  case  undershoots  the  peak  at  X = 1 by  about  20%  which  is  about 
1.0  dB.  The  periodic  estimate  is  much  improved  in  terms  of  aliasing  for  this  sample  rate  but 
is  still  poor  for  frequencies  between  X = 1.0  and  X = 1.5.  The  estimate  is  biased  in  amplitude 
for  most  of  the  X but  for  portions  of  the  spectral  density  the  estimate  is  reasonable.  If  we 
had  extended  the  frequency  axis  of  the  graphs  for  the  periodic  case  another  pair  of  peaks 
would  erroneously  appear  every  multiple  of  27rp.  This  would  not  occur  for  the  Poisson  esti- 
mate. The  third  sampling  rate,  p = ^ = 3.2  rad/s,  was  used  to  generate  the  estimates  shown  i 

in  Figure  9.  In  this  case,  for  X < 1 .5,  the  periodic  estimate  is  very  good  but  for  X > 1 .5  the 
aliasing  error  dominates.  The  Poisson  estimate  gives  a good  representation  of  the  true  spec- 
tral density  over  the  entire  range  of  the  graph. 

C.  Two  Closely  Spaced  Signals 

This  case,  as  described  earlier,  was  generated  by  adding  the  time  series  of  two  narrow- 
band  signals  with  different  center  frequencies.  The  center  frequency  spacing  was  determined 
in  such  a way  that  the  null  between  the  two  peaks  is  located  at  the  common  half  power 
point  (-3  dB)  of  the  two  peaks.  The  values  of  the  parameters  in  equation  (III-28)  used  for 
the  simulation  are  the  same  as  in  the  narrowband  case  with  the  addition  of  the  parameter 
Xj  = 1.334.  Because  the  bandwidth  (III-9)  is  dependent  on  the  center  frequency  the  band- 
width and  the  height  of  the  two  peaks  differ  slightly.  Two  sampling  rates  were  used  in  eval- 
uating the  estimators  for  this  spectral  density; 

p = 0=2ira  rad/s  (IV-3a) 

p = /3  = 3.2  rad/s.  (IV-3b) 

These  two  sampling  rates  provide  an  example  of  sampling  well  below  the  two  peaks  and 
another  of  sampling  above  them  in  frequency.  Other  parameters  used  for  tliis  spectral  den- 
sity are  identical  to  the  narrowband  case  (i.e.  N = 1 000,  a = 2. 1 /22).  The  results  for  this 
spectral  density  are  presented  in  Figures  10  and  1 1.  Figure  10  shows  the  results  for  the 
sample  rate  (27ra  rad/s).  The  randomly  sampled  estimate  (top  of  Figure  10)  shows  fairly 
I large  oscillations  with  frequency  but  does  succeed  in  providing  a fair  estimate  of  the  spectral 

; density.  It  is  of  significance  that  even  though  the  crossover  point  between  the  two  peaks  is 
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Figure  8a.  Poisson  spectral  estimation  of  “narrowband”  spectral  density  sampled  at 
P = 2.2  rad/s  with  Mj^  = 55. 


Figure  8b.  Periodic  spectral  estimation  of  “narrowband”  spectral  density  sampled  at 
p = 2.2  rad/s  with  Mj^j  = 50 
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Figure  10a.  Poisson  spectral  estimation  of  two  nearby  “narrowband”  signals  sampled  at 
P=  Ina  ~'0.6  rad/s  with  Mj^j  = 55. 
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. Figure  10b.  Periodic  spectral  estimation  of  two  nearby  “narrowband”  signals  sampled  at 

p = 2»t  a ~0.6  rad/s  with  M]y  = 20- 
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Figure  I la.  Poisson  spectral  estimation  of  two  nearby  “narrowband  signals  sampled  at 
0 = 3.2  rad/s  with  = 60- 


FREQUENCY  (rad/«) 

Figure  I lb.  Periodic  spectral  estimation  of  two  nearby  “narrowband  signals  sampled  at 
p = 3.2  rad/s  with  Mjyj  = 55. 
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only  3 dB  the  estimate  does  indeed  separate  the  peaks.  For  the  periodically  sampled  estimate  j 

there  is  no  resemblance  of  the  estimate  to  the  true  spectral  density.  In  this  case  the  aliased 

spectral  density  is  not  even  well  estimated.  It  should  be  noted  that  the  amplitude  axis  has 

been  scaled  down  in  the  graphs  to  accommodate  the  large  aliasing  effects  from  the  periodic 

estimates.  Figure  1 1 , which  is  for  a sample  rate  of  3.2  rad/s,  is  improved  for  both  the 

random  case  and  the  periodic  case.  The  random  estimate  demonstrates  less  variance  but 

otherwise  has  changed  very  little.  This  shows,  as  in  the  other  spectral  densities  simulated,  j 

that  the  Poisson  technique  does  not  depend  strongly  on  the  sampling  rate.  On  the  other  i 

hand,  the  periodic  estimate  shown  in  Figure  10 is  much  improved,  particularly  for  X < 1 .5.  | 

Still,  for  the  estimate  at  higher  frequencies  in  Figure  1 1 the  aliasing  error  drastically  dis-  ; 

torts  the  estimate.  i 

D.  Data  Integrity  | 

In  a simulation  study  where  synthetically  generated  data  are  used  the  question  concern- 
ing the  quality  of  the  data  arises.  The  theory  assumes  the  Poisson  sample  spacings  (I1I-6)  and 
the  data  X(t)  are  uncorrelated  from  sample  to  sample.  In  addition,  the  Poisson  sample  spac-  i 

ings  and  the  data  are  assumed  not  to  be  cross  correlated.  Unfortunately,  for  a simulation  of 
this  type  that  is  not  possible. 

A random  number  generator*  was  used  to  generate  exponentially  distributed 
(mean  = 1)  random  numbers  lor  the  Poisson  sample  times  {t)^  }|^=)  uid  normal 

(0,  1)  random  numbers  which  were  used  as  a basis  for  simulating  the  data  . Sample  | 

auto-  and  cross-correlations  for  the  random  numbers  used  in  the  simulation  were  calculated.  ; 

I 

The  sample  correlations  resulting  from  random  numbers  used  in  the  first  order  Clauss-Markov 
spectral  density  for  lags  to  ±10  were  calculated  and  the  more  significant  values  are  given  in  | 

Table  1.  For  the  autocorrelations,  normalized  to  1 at  zero  lag,  the  correlations  were  consid-  | 

ered  acceptable  with  the  exception  of  {a|^}  for  lag  ± 1 . Similarly  the  crosscorrelation  values  ] 

under  1 .0%  were  considered  good.  It  is  possible  that  numbers  with  smaller  sample  correlation  j 

values  could  have  been  generated  with  successive  computer  runs,  but  numbers  with  sample 
correlations  of  this  order  were  considered  to  be  representative  of  the  random  number  gener- 
ator and  therefore  the  results  of  this  simulation  could  easily  be  duplicated. 

Another  approach  would  have  been  to  generate  many  sets  of  numbers  and  compute 
spectral  estimates  for  each  set  of  numbers  and  average  them  together.  This  would  have  been 
very  costly  and  since  this  simulation  compares  two  different  spectral  estimation  techniques, 
this  was  not  necessary. 

Sample  correlation  values  for  random  numbers  used  in  the  other  spectral  densities 
simulated  were  very  similar  to  those  in  Table  1 and  are  not  presented  here. 

*A  software  random  number  generator  from  the  International  Mathematics  and  Statistics  Library  ( IMSL) 
was  used  for  the  random  numbers.  The  UNIVAC  1 1 10  computer  was  used  in  the  simulation. 
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V.  CONCLUSIONS 

The  objective  of  this  work  was  to  provide  a fair  comparison  between  two  techniques 
for  estimating  the  spectral  density  of  a continuous  time  stationary  process. 

The  main  feature  of  the  Poisson  spectral  estimation  technique  is  its  ability  to  estimate 
the  continuous  time  spectral  density  0{X)  rather  than  the  aliased  discrete  spectral  density 

It  is  this  property  that  eliminates  the  problem  of  spectral  aliasing  which  is  inherent  in 
some  degree  in  all  real  world  spectral  analysis. 

From  the  results  presented  in  this  paper  it  can  be  concluded  that  the  Poisson  tech- 
nique out-performs  the  periodic  technique  in  a wide  variety  of  circumstances.  This  is  signifi- 
cant because  the  periodic  technique  is  the  standard  method  used  in  spectral  analysis.  In  the 
first  order  Gauss-Markov  case  the  Poisson  technique  was  shown  to  yield  good  estimates  of  a 
“wide  band”  low  pass  spectral  density.  The  successful  estimation  of  the  more  complicated 
“narrowband”  spectral  density  demonstrated  the  Poisson  technique’s  ability  to  handle  spec- 
tral estimates  of  the  “narrowband”  class.  Finally,  in  the  last  case  considered,  the  Poisson 
technique  showed  the  ability  to  distinguish  two  “narrowband”  signals  closely  spaced  in  fre- 
quency. These  three  spectral  densities  were  chosen  for  the  simulation  largely  because  they 
represent  the  type  of  spectra  typically  encountered  in  ocean  data. 

The  comparison  of  the  Poisson  and  periodic  techniques  for  the  spectral  densities 
studied  showed  the  Poisson  to  clearly  excel  in  the  following  cases: 

a.  When  the  spectral  density  is  not  bandlimited. 

b.  When  the  shape  of  the  spectral  density  is  not  known  a priori. 

c.  When  the  sampling  rate  is  constrained  to  be  below  that  of  Nyquist. 

If  the  signal  to  noise  ratio  and  noise  power  are  known  beforehand  the  optimal  sampling  rate 
for  them  may  be  computed  and  the  Poisson  technique  used  to  provide  spectral  estimates 

The  main  disadvantage  of  the  Poisson  technique  is  that  it  is  more  difficult  to  imple- 
ment than  the  periodic  technique  since  FFT  (fast  Fourier  transform)  type  procedures  cannot 
be  used.  Development  of  techniques  to  efficiently  extract  Poisson  distributed  data  samples 
from  continuous  time  data  and  compactly  perform  the  spectral  estimates  would  make  the 
Poisson  technique  more  appealing  for  practical  applications. 

The  Poisson  theory  applies  to  other  areas  of  signal  processing  as  well  as  spectral  den- 
sity estimation.  This  could  result  in  using  the  Poisson  technique  in  applications  such  as  cross 
spectral  estimation,  estimation  of  the  coherence  function,  and  antenna  array  beamforming. 
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APPENDIX 

In  general,  if  X has  a rational  spectral  density  function,  i.e. 


Q(iX) 


where 


a)  P(\)  = ^ bj  (iX)* 

n 

b)  Q(X)  = 2 aj  (iX)j 

j=0 

and  a,^  = I , with  m < n,  then  X(t)  formally  satisfies  the  following  stochastic  differential 
equation: 

X(")(t)  + an_i  x('^-‘>(t)  + . . . + aoX(t)  = (t)  + . ..  + boU(t).  (A-2) 

U(t)  is  a stationary  white  noise  process  with  E U(t)  U*(s)  = 6(t-s). 

' Now,  (A-1 ) can  be  represented  as 


X(t)  = FX(t)  + GU(t) 


where 


X(t)  = 


X,(t) 

XtU) 


l_X,^(t)  J 


with 


X,(t)  = X(t) 

X2(t)  = X,(t)  + a„_,  X(t)-b„_,  U(t) 


Xn(t)  = X„(t)  + a,  X(t)-b,  U(t) 
The  solution  of  (A-l ) is  given  by 

t 


X(t)  = <t(t,  to)X(to)+  J 4>(t,T)G(r)U(T)dT 


t^ 


(A-4) 


where  the  state  transition  matrix  <t>(t,tQ)  satisfies  the  (deterministic)  differential  equation 
^4>(t,r)  = F^(t,T) 
whose  solution  is  given  by 

4>(t,T)  = (A-5) 

Equation  (A-4)  is  known  as  the  dynamic  representation  of  the  random  process  X. 


